%% experimental data
load('/Users/nzhao/Code/diffuse-spin/APP/PumpTrans_sendtoDW.mat')
Pin = cell2mat(inputPower);
Pout = cell2mat(outputPower);

%% param
param = struct( 'cellSize',          1.0, ...     % cm
                'cellTemperature',   83, ...   % degree Celsius
                'magneticField',     21.5e-6, ... % Tesla
                ...
                'Rb85Abundance',     0.7217, ...  % natural abundance 0.7217
                'Rb87Abundance',     0.2783, ...  % natural abundance 0.2783
                'N2Pressure',        450, ...     % Torr
                'Xe129Pressure',     14, ...       % Torr
                'Xe131Pressure',     4, ...      % Torr
                ...
                'PumpPower',         100e-3, ...   % W
                'PumpPolarization',  [1, 1j], ... % circularly polarized
                'PumpDirection',     [0 0 1], ... % z direction
                'PumpWaist',         0.25, ...    % cm
                'PumpDetuning',      -4.5e3, ...  % MHz
                ...
                'ProbePower',        1e-3, ...    % W
                'ProbePolarization', [1, 0], ...  % linearly polarized
                'ProbeDirection',    [1 0 0], ... % x direction
                'ProbeWaist',        0.1, ...     % cm
                'ProbeDetuning',     -20.0e3...   % MHz
                );

%% initialization
nnmrg = Table0.NumericalNMRG();
nnmrg.initializeSystem(param);
%%
pumpAbsorption = nnmrg.GetCellAbsorption('pump', 'frequencyRange', -30e3:30e3, 'isPlot', true);
probeAbsorption = nnmrg.GetCellAbsorption('probe', 'frequencyRange', -30e3:30e3, 'isPlot', true);

%%
rateRelax = 5e3;
zList = 0:0.01:param.cellSize; xList = 0:0.01:param.cellSize/2; [xx, zz] = meshgrid(xList, zList);
res=cell2mat(arrayfun(@(x) nnmrg.GetBeamPropagation(x, 0, zList, rateRelax, 'output', 'rate'), xList, 'UniformOutput', false));

f2 = figure; ax2 = axes('Parent', f2);
contourf(ax2, xx, zz, res/2/pi/1e3, 20); colorbar(ax2);

%%
temperatureList = [83.8, 100.2, 106.6, 111.3, 115.1, 118.1, 120.9];
power_in = [0.01e-3, 1e-3:2e-3:9e-3, 10e-3:10e-3:180e-3]';

[transmission, transmission1] = arrayfun(@(t) compare_transmission(nnmrg, param, rateRelax, power_in, t), temperatureList, 'UniformOutput', false);
%%
f=figure; ax=axes('Parent', f);
cla(ax); hold(ax, 'on');
plot(ax, power_in*1e3, cell2mat(transmission), '-', 'linewidth', 2); 
plot(ax, Pin, Pout./Pin/0.72, 'o'); 
grid(ax, 'on'); hold(ax, 'off'); 
xlim(ax, minmax(power_in'*1e3)); legend(ax, TempLegendcell, 'Location', 'best');

%%
function [transmission, transmission1] = compare_transmission(nnmrg, param, rateRelax, power_in, temperature)
    param.cellTemperature = temperature;
    transmission = arrayfun(@(p) calc_transmission(nnmrg, param, rateRelax, p), power_in(:));
    transmission1 = arrayfun(@(p) calc_transmission1(nnmrg, param, rateRelax, p), power_in(:));
end
%%
function res = calc_transmission(nnmrg, param, rateRelax, power_in)
    param.PumpPower = power_in;
    nnmrg.initializeSystem(param);
    power_in1 = integral2(@(x, y) nnmrg.GetBeamPropagation(x, y, 0.0, rateRelax, 'output', 'intensity'), -.75, .75, -.75, .75);
    power_out = integral2(@(x, y) nnmrg.GetBeamPropagation(x, y, param.cellSize, rateRelax, 'output', 'intensity'), -.75, .75, -.75, .75);
    res = power_out/power_in1;
end

function res = calc_transmission1(nnmrg, param, rateRelax, power_in)
    param.PumpPower = power_in;
    nnmrg.initializeSystem(param);
    intensity_in = nnmrg.GetBeamPropagation(0, 0, 0, rateRelax, 'output', 'intensity');
    intensity_out = nnmrg.GetBeamPropagation(0, 0, param.cellSize, rateRelax, 'output', 'intensity');
    res = intensity_out/intensity_in;
end
    

